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Description of the Chemical Reaction Path in the HCO Molecule: 
A Combined Configuration Interaction and Tight-Binding Approach 
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It is demonstrated that the reaction path for a polyatomic molecule (applied to the HCO molecule) 
is easily calculated via geometry-independent tight binding Hamiltonian fitted to accurate ab-initio 
configuration interaction (CI) total energies. This Hamiltonian not only reproduces the CI calcula- 
tions accurately and efficiently, but also effectively corrects any CI energies happening to erroneously 
converge to excited states. 

PACS numbers: 31.10.+Z 31.15.-p 31.50.-x 82.20.Kh 
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The question. The determination of the reaction path 
in a chemical reaction needs the detailed knowledge of 
the pertinent potential energy surface (PES) (diabatic 
or adiabatic). This is a formidable task because (i) the 
PES is a multi-dimensional surface, impossible to be ab- 
initio calculated at every point in the degrees of free- 
dom (DOF) space, and interpolation is necessary (the 
most accurate known detailed multi-dimensional PES is 
that of H3 interpolated from 71969 ab-initio DOF points 
0). (ii) Because the standard ab-initio calculations of 
the many electron problem in the Born-Oppenheimer ap- 
proximation, based on the variational principle, (accu- 
rate ab-initio configuration interaction (CI) calculations) 
yield only adiabatic curves, and, more importantly, be- 
ing iterative, sometimes converge to undesirable states 
also c.f. below). Yet, such calculations may be in- 
hibitively time consuming. For the ground state, the time 
problem is already traditionally overcome via the density 
functional theory (DFT) @, which self-consistently ap- 
proximates the many electron by a one-electron problem. 
However, DFT calculations sometimes fail to explain ex- 
perimentally observed features of the PES [4( . Thus, the 
accurate CI calculations are more or less indispensable, 
even if performed in a rather limited, but representa- 
tive, set of molecular geometries. Therefore, a reliable 
interpolation scheme for the pertinent PES, based on CI 
calculations, and overcoming the problem of wrong CI 
convergence, is desirable. 



The purpose. It is shown that such an interpo- 
lation scheme is possible, based on a spin-polarized 
H geometry-independent @ Slater - Koster (SK) 
parametrization pj of ab-initio CI total energies |jg. As 
a demonstration, the method is applied to a triatomic 
molecule of chemical kinetics interest, HCO, which is an 
intermediate radical in the generation of a primary ion 



during hydrocarbon combustion: 
0( 3 P) + CH(X 2 U and/or a 4 XT) 



[hco]*Ca') 

HCO+ + e~ 



The reaction of 0( 3 P) with CH(X 2 n, a 4 ST) is known 
experimentally Q to generate the HCO + cation via au- 
toionization of some state (or states) of the intermedi- 
ate HCO radical upon interaction with some vibrational 
level of the ion. The first step toward computations of 
such interactions is the construction of the potential en- 
ergy surface (PES) of the states with low (or no) barrier, 
through which a reaction at the experimental tempera- 
ture can proceed. Such a state (without a barrier) is the 
HCO(X 2 A') state ^(j > used here to test the interpolation 
scheme on a molecular system. The reaction path of the 
formation of HCO (in 2 A' symmetry) is also computed 
using the interpolated PES. 

The procedure. First several (724 - compared to 71969 
of H 3 1]) accurate CI total energies, based on (less 
accurate) multi-configuration self-consistent field (MC- 
SCF) orbitals, are calculated at selected geometries of 
the H,C,0 atoms in the A' symmetry of the Cs group. 
Most of them (508) are fitted to the interpolation scheme, 
the remaining serving to check the quality of the fit. 
For the fit a non-orthogonal spin-polarized tight bind- 
ing (TB) Hamiltonian is formed, whose matrix elements, 
along with those of the overlap matrix, are expressed 
as functions of the bond direction, according to the SK 
scheme Q, and of the bond length, according to the 
Naval Research Laboratory (NRL) technique [||, he.: 
The functions are generally polynomials of the inter- 
atomic distance, within exponential envelopes, the co- 
efficients and the exponents being varied as parameters. 
For two adiabatic states near some (avoided) crossing the 
TB Hamiltonian naturally produces two diabatic PESs 
in nearby extrapolation, and predicts to which diabatic 
PES, ground-state or excited, nearby CI energies belong. 
Among these, the appropriate ones can be used to ex- 
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tend the fit beyond the (avoided) crossings, around which 
two sets of parameters are needed for the two PES's. If 
it happens, as with HCO, that the ground and excited 
state energies beyond the crossing lie close to each other, 
the adiabatic PES can be fitted as well, with comparable 
accuracy. 

Finally, by using at each point of the DOF space the 
lowest lying TB-fitted PES, the adiabatic path can be 
found: For each value of a desired degree of freedom (in 
our case for each C-0 distance) the energy minimum is 
searched in the space of the remaining degrees of 
freedom (C-H distance and H-C-0 angle). Having the 
parametrized tight binding Hamiltonian, any property 
can be trivially computed. 

Methodology. For the CI energies the correlation con- 
sistent aug-cc-pVTZ basis set was used ^2 El m con ~ 
junction with the complete active space self-consistent 
field (CASSCF) + 1 + 2 multi-reference CI method 
(MRCI) employed in the MOLPRO package @ (the four 
electrons in the Is orbitals of C and O were unexcited). 
The CASSCF calculations were state-averaged, and the 
active space was limited to the 9 valence orbitals among 
which the remaining 11 electrons were distributed. In 
the subsequent MRCI calculations the uncontracted con- 
figurations were around 50 million internally contracted 
to about one million. Calculations between C-0 dis- 
tances of 1.7 and 6 bohr were done for several H-C- 
O angles between 50° and 180° and several C-H dis- 
tances between 1.7 and 4.5 bohr, most around the C- 
H equilibrium distance of 2.12 bohr. The three lowest 
roots of the secular equation were computed to increase 
the accuracy of the calculation. By an analytic gradi- 
ent optimization at the MCSCF level, an approximate 
(MCSCF) equilibrium geometry was found at the DOF 
space point \r HC , r CO , 9 B -c-o) = (2.12, 2.2, 126°) (in 
a.u.). Because it is not evident whether the aforemen- 
tioned points are beyond any avoided crossing, where 
the role of the ground and the excited states would be 
interchanged, first several DOF points near equilibrium 
were obtained by employing a generalization of the 3- 
dimensional sphere to the generally multi-dimensional (in 
this case also 3-dimensional) DOF space: X\ = rj/fj — 1, i 
= {HC, CO}, £3 = 9/8— 1, where generally for n degrees 
of freedom, points belonging to a ri-dimensional hyper- 
sphere of radius r and center (a^, i = l,...,n) are obtained 
by 

x n -i - x n -i = r sinO n cos9 n -\ (1) 

X\—X\ — r sin8 n sin9 n -\...cos6i 

where the 1st B\ = or 180°, the two points of a "1- 
dimensional sphere" , and the other < 9i < 180° are the 
"azimuthal" hypersphere angles (incidentally, a variable 
dimensional do-loop code was invented, needed to treat 
any larger molecule). Thus, first points with small r were 
fitted, and gradually the fit was extended to more remote 



DOF points. 

The formalism of the NRL geometry - independent TB 
parametrization is described in detail in Ref . ; here an 
essential summary is only presented. The total energy is 
written as 

E[n(f)} = /(^^) e * • + F Hr)} 

i ; s=l,2 

- £ k^t^)^* (2) 

i ; s=l,2 

where fix) = 1/(1 + e x ), T=0.005 mRy, and 

ei' s =e is +V ; fx' = fX + V ; V = F[n(r)]/N e (3) 

with N e = J^i ■ s=i.2 /((A* _ e i s)/T) being the number 
of electrons, i counts the states, s = 1, 2 counts the spin. 
Since the total energy is independent of the choice of 
zero of the potential, the shift Vq is sufficient to be deter- 
mined by the requirement that e/ s are the eigenvalues of 
the generalized eigenvalue problem (H — S e{ s ) ipi s = 0, 
where H is the TB Hamiltonian and S is overlap matrix 
in an atomic s- and p-orbital basis representation {<fi a }- 
Thus, a non-orthogonal TB calculation uses on-site, hop- 
ping and overlap parameters. Demanding that only the 
on-site SK parameters are affected by the shift Vq, for 
atom / in a spin-polarized structure the matrix elements 
are expressed as 

h\ s = jzH ns £T ; l = s,P (4) 

71=0 

where 

W . = XV AJ "«" /(^LZi?.) (5 ) 

is a generalized pair potential ("density"), with Rq = 15 
bohr, r c = 0.5 bohr, Ri j is the internuclear distance be- 
tween atoms / and J, I(J) denote the type of atom on 
the site I (J) while A/ j , depending on the atom type, 
and bf n s are the on-site NRL geometry-independent pa- 
rameters (GIP). It is found sufficient to keep hopping and 
overlap parameters spin independent, of the form 

2 

P,(R) = (£c in R n ) e~Z R f{— (6) 

where 7 indicates the type of interaction (i.e. ssa, spcr, 
ppo~, ppn and psa). The NRL GIPs are c 7 „ and <? 7 , R 
is the interatomic distance, and Rq and r c are as in eq. 

_ 

Within the context of the NRL code [6j , written pri- 
marily for solids, the molecule was treated as a base to a 
large cubic lattice unit cell (lattice constant = 100 a.u.) 
ensuring vanishing interaction between atoms in neigh- 
boring cells. Thus, the PES was described in terms of 
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the following NRL GIPs for each spin polarization. On- 
site: s: H, C, O, (H depending on C), (C on H), (H on 
O), (O on H), (C on O), and (O on C); p: C, O, (C on 
H), (O on H), (C on O), and (O on C). Hopping and 
overlap parameters: ssct: H-C, H-O, C-O; spc: H-C, H- 
O, C-0 and O-C (denoted as pscr); ppcr and pp7r: C-O. 
For HCO, since similar atoms are well separated, the H- 
H, C-C and 0-0 parameters vanish. We fitted 508 CI 
points and checked the resulting PES against 216 more 
CI energies not included in the fit. The error was less 
than 10~ 3 a.u., which is within the ab-initio PES accu- 
racy (starting from different initial guesses the MCSCF 
calculation may converge to slightly different results by 
1CP 3 a.u.). To ensure obtaining physically meaningful 
TB parameters, for a very limited number of molecular 
geometries the Hamiltonian eigenvalues were also fitted, 
while the total energy was fitted for all 508 structures. 

Finally, for the reaction path we used a non-linear en- 
ergy minimization technique employing Powell's conju- 
gate directions method |l5J modified to be restricted to 
closed intervals of the DOF space • 

For comparison, each of the 724 ab-initio CI calcu- 
lations needs 3 hours of CPU time, each n-dimensional 
hypersphere radius r-increase, to fit more remote points 
(with 10 such hypersphere radial extensions all points 
can be covered) needs 2-3 hours, and each 2-dimensional 
energy minimization, using the final TB parameters (i.e. 
the reaction path determination), needs a few seconds. 

Results. The fitted TB Hamiltonian could predict cor- 



rectly total energy curves for points not included in the 
fit as shown for example in Fig. ^ Since it produces 
naturally the diabatically extended branch of the en- 
ergy, it could distinguish to which adiabatic state near an 
avoided crossing the CI values belong. Classifying such 
CI points may sometimes be misleading or unrecogniz- 
able by mere observation of the MCSCF orbitals. An ex- 
ample is shown in Fig. [21 However, the most impressive 
aspect was that we realized, through the fit, that at some 
points (about 10 in 700) the CI calculation had converged 
to excited energies (which ought to be disregarded, oth- 
erwise they would destroy the fit). An example is given 
in Fig. [3J Finally, Fig. 0] shows the reaction path for the 
formation of HCO, as HC approaches O. For a triatomic 
molecule the figure contains the whole information: For 
each C-O distance the minimum energy and the corre- 
sponding C-H distance and H-C-0 angle are displayed. 
At large C-O distances, O is more attracted toward H, 
but, in approaching equilibrium, O binds mainly with C, 
the H-C-0 angle gradually becoming ~ 122° (represent- 
ing the CI value). Around equilibrium (c.f. Table QJ, 
the angle changes slightly monotonically by 1-2°, but be- 
cause, in increasing the C-O distance, the C-H distance 
decreases, predominantly an antisymmetric stretching vi- 
bration occurs. To our knowledge there is no experimen- 
tal confirmation of the reaction path of this intermediate 
molecule. 
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Fig 1a 
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FIG. 1: Predicted total energy E in a.u. (Above:) vs C-O 
distance for C-H distance = 3.01 bohr, and various H-C-O 
angles. (Below:) vs C-H distance for various C-O distances, 
and H-C-O angle = 100°. 
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Fig.2 

C-H= 1.71 a.u 
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FIG. 2: The CI point A (excited) in E vs C-O distance for C- 
H distance = 1.71 bohr and H-C-O angle =180° is predicted 
by the fit to belong to the diabatic branch of the curve beyond 
the avoided crossing. (Inclusion of the lower value to the fit, 
destroys it.) 




C-0 distance (1 .80 to 4.02 a.u.) 



FIG. 3: The CI points A and B clearly belong to the excited 
state as shown by the TB prediction. The CI calculation could 
not converge to the correct values. The discontinuity can be 
verified by observing the corresponding MCSCF orbitals. 
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FIG. 4: The reaction path for the formation of HCO. Details 
are described in the text. 



